//============================================================================
//                                  I B E X
// File        : ibex_LSmear.cpp
// Author      : Ignacio Araya, Bertrand Neveu
// License     : See the LICENSE file
// Created     : May, 25 2017
// Last Update : May, 3 2019
//============================================================================

#include "ibex_LSmear.h"
#include "ibex_ExtendedSystem.h"
#include "ibex_Random.h"
#include "ibex_OptimLargestFirst.h"

using std::pair;

namespace ibex {

LSmear::~LSmear() {
	delete mylinearsolver;
}


  LPSolver::Status LSmear::getdual(IntervalMatrix& J, const IntervalVector& box, Vector& dual) const {

    int  _goal_var = goal_var();
    bool minimize=true;
    if (_goal_var == -1){
      _goal_var = RNG::rand()%box.size();
      minimize=RNG::rand()%2;
    }

	// The linear system is created
	mylinearsolver->clear_constraints();
	mylinearsolver->set_bounds(box);
	mylinearsolver->set_bounds(_goal_var, Interval(-1e10,1e10));

	int* nb_lctrs = new int[sys.f_ctrs.image_dim()]; /* number of linear constraints generated by nonlinear constraint*/

	for (int i=0; i<sys.f_ctrs.image_dim(); i++) {

		if (J[i].is_unbounded()) {
			nb_lctrs[i]=0;
			continue;
		}

		Vector row1(sys.nb_var);
		Interval ev(0.0);
		for (int j=0; j<sys.nb_var; j++) {
			row1[j] = J[i][j].mid();
			ev -= Interval(row1[j])*box[j].mid();
		}
		ev+= sys.f_ctrs.eval(i,box.mid()).mid();

		nb_lctrs[i]=1;
		if (i!=goal_ctr()) {
		  if (sys.ops[i] == LEQ || sys.ops[i] == LT) {
		    mylinearsolver->add_constraint( row1, sys.ops[i], (-ev).ub());
		  }
		  else if (sys.ops[i] == GEQ || sys.ops[i] == GT)
				mylinearsolver->add_constraint( row1, sys.ops[i], (-ev).lb());
		  else { //op=EQ
				mylinearsolver->add_constraint( row1, LT, (-ev).ub());
				mylinearsolver->add_constraint( row1, GT, (-ev).lb());
				nb_lctrs[i]=2;
		  }
		}
		else if (goal_to_consider(J,i) ||sys.f_ctrs.image_dim()==1 )
		  mylinearsolver->add_constraint( row1, LEQ, (-ev).ub());
		else // the goal is equal to a variable : the goal constraint is useless.
		  nb_lctrs[i]=0;

	}


	//the linear system is solved
	LPSolver::Status stat=LPSolver::Status::Unknown;

	mylinearsolver->set_cost(_goal_var, (minimize)? 1.0:-1.0);

	stat = mylinearsolver->minimize();

	if (stat == LPSolver::Status::Optimal) {
		// the dual solution : used to compute the bound
		Vector dual_solution(mylinearsolver->nb_rows());  //
		dual_solution = mylinearsolver->not_proved_dual_sol();

		int k=0; //number of multipliers != 0
		int ii=0;
		// writing the dual solution in the dual Vector of dimension +
		for (int i=0; i< sys.nb_var; i++)
			dual[i]=dual_solution[i];
		for (int i=0; i<sys.f_ctrs.image_dim(); i++) {
			if (nb_lctrs[i]==2) {
				dual[sys.nb_var+i]=dual_solution[sys.nb_var+ii]+dual_solution[sys.nb_var+ii+1]; ii+=2;
			}
			else  if (nb_lctrs[i]==1) {
				dual[sys.nb_var+i]=dual_solution[sys.nb_var+ii]; ii++;
			}
			else {
				dual[sys.nb_var+i]=0.0;
			}
			if (std::abs(dual[sys.nb_var+i])>1e-10) k++;
		}

		if(k<2) { stat = LPSolver::Status::Unknown; }
	}
	delete[] nb_lctrs;
	return stat;
}


int LSmear::var_to_bisect(IntervalMatrix& J, const IntervalVector& box) const {
  int lvar = -1;

	if (box.is_unbounded()) {
		return SmearSumRelative::var_to_bisect(J, box);
	}
	//Linearization
	LPSolver::Status stat = LPSolver::Status::Unknown;

	Vector dual_solution(sys.nb_var+sys.f_ctrs.image_dim());

	if (lsmode==LSMEAR_MG) { //compute the Jacobian in the midpoint
		IntervalMatrix J2(sys.f_ctrs.image_dim(), sys.nb_var);
		IntervalVector box2(IntervalVector(box.mid()).inflate(1e-8));
		//		IntervalVector box2(IntervalVector(box.random()));
		box2 &= box;

		sys.f_ctrs.jacobian(box2,J2);
		stat = getdual(J2, box, dual_solution);

	} else if (lsmode==LSMEAR) {
		stat = getdual(J, box, dual_solution);
	}


	if (stat == LPSolver::Status::Optimal) {
		double max_Lmagn = 0.0;
		int k=0;

		for (int j=0; j<sys.nb_var; j++) {
			Interval lsmear=Interval(0.0);
			if ((!too_small(box,j)) && (box[j].mag() <1 ||  box[j].diam()/ box[j].mag() >= prec(j))){
				lsmear=dual_solution[j];

				for (int i=0; i<sys.f_ctrs.image_dim(); i++){
				  lsmear += dual_solution[sys.nb_var+i] * J[i][j];
				}
			}

			lsmear*=(box[j].diam());

			if (lsmear.mag() > 1e-10  && (j!=goal_var() || mylinearsolver->minimum().mid() > box[goal_var()].lb() )) {
				k++;
				if (lsmear.mag() > max_Lmagn) {
					max_Lmagn = lsmear.mag();
					lvar = j;
				}
			}
		}

		if (k==1 && lvar==goal_var()) { lvar=-1; }
	}
	if (lvar==-1) {
	  //	  std::cout << "ssr " << std::endl;
	  lvar=SmearSumRelative::var_to_bisect(J, box);
	}
	//	std::cout << "lsmear " << lvar << std::endl;
	return lvar;
}


} /* namespace ibex */
